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We propose a general method of using the Fokker-Planck equation (FPE) to link the Monte- 
Carlo (MC) and the Langevin micromagnetic schemes. We derive the drift and diffusion FPE 
terms corresponding to the MC method and show that it is analytically equivalent to the stochastic 
Landau-Lifshitz- Gilbert (LLC) equation of Langevin-based micromagnetics. Subsequent results 
such as the time quantification factor for the Metropolis MC method can be rigorously derived from 
this mapping equivalence. The validity of the mapping is shown by the close numerical convergence 
between the MC method and the LLG equation for the case of a single magnetic particle as well 
as interacting arrays of particles. We also found that our Metropolis MC is accurate for a large 
range of damping factors a, unlike previous time-quantified MC methods which break down at low 
a, where precessional motion dominates. 
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With the rapid advance of computing resources, Monte 
Carlo (MC) methods have become a powerful tool in 
many fields ranging from the physical sciences to fi- 
nance and sociology 0, 0, 13 • The flexibility of Monte 
Carlo is due to its abstract formalism which can be real- 
ized in an almost infinite number of ways. Increasingly, 
MC methods are being implemented in the stochastic 
micromagnetic modeling of magnetic nanostructures . 
Stochastic micromagnetic modeling has important prac- 
tical implications, e.g. in predicting the storage lifetime 
of hard-disk magnetic media ^ Traditionally, the 
dynamics of magnetic moments are also modeled in the 
Langevin scheme, using the stochastic Landau-Lifshitz- 
Gilbert (LLG) equation "T]. Langevin-based micromag- 
netics constitutes a formidable computational method, 
because of its ease of use and close correspondence to ac- 
tual experimental data in previous literatures 1^. How- 
ever, it has certain limitations which may be overcome by 
MC methods. For instance, MC methods can accommo- 
date the long-time magnetization relaxatio n dy namics of 
large-scale arrays of magnetic grains |3, llfl [HI i which is 
practically unfeasible to be modeled using the stochas- 
tic LLG equation. On the other hand, MC schemes 
have the drawback of having its time calibrated in MC 
steps, instead of physical time units. Thus, both MC and 
the Langevin approach are useful computational methods 
in micromagnetics, with complementary strengths and 
drawbacks. Hence, it is very important to devise a gen- 
eral way of mapping one method to the other, and vice 
versa. 

Early efforts to link MC to LLG were done by Nowak et 
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al. They focused on deriving a time quantification 

factor to relate one MC step (MCS) to real physical time 
unit used in the LLG equation. Recently, we also pro- 
posed another time-quantifiable MC (TQMC) method 
which involves the determination of macroscopic density 
of states, and the use of the Master equation for time evo- 
lution. This method is applicable in simulating extremely 
long time magnetization reversal process 13]. The effect 
of precession on Nowak's time-quantification was inves- 
tigated by Chubykalo et al. 14]. They concluded that 
Nowak's time quantification of MC breaks down in the 
low damping case, in presence of an oblique external field, 
due to the influence of (athermal) precessional motion. 

In this paper, our approach is more thorough and fun- 
damentally different from all previous works. We propose 
a systematic way of using the Fokker-Planck equation 
(FPE) to map MC to LLG dynamical equation. The 
physical background of using FPE to describe stochas- 
tic dynamics has been well established ^-S- i^i the 
case of a particle under the influence of a one-dimensional 
potential T^. The outline of our scheme is as follows. 
We consider a single isolated particle and then general- 
ize to an interacting array of particles. First, we develop 
a MC method that has the stochastic dynamics of the 
LLG equation. This method is a hybrid Metropolis MC 
scheme which combines the random displacement of spins 
about a cone and a suitably-sized precessional step. The 
random displacement models the thermal fluctuation and 
the precessional step accounts for the precessional term 
in the LLG equation. From the drift-diffusion picture of 
the general Fokker-Planck equation, we can then obtain 
the FP coefficients of this hybrid Metropolis MC method, 
and compare them with previously obtained FP coeffi- 
cients of the stochastic LLG equation. The comparison 
shows an exact equivalence of the two EPFs in the limit 
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of small MC step size. From this comparison, we obtain 
the time quantification factor to relate one MC step to 
real time unit and the required step size for the preces- 
sion, which allows us to map the MC results to that of 
the LLC equation. Generalization to an array of inter- 
acting particles will be shown numerically in the latter 
part of this paper. 

We consider an isolated single domain magnetic par- 
ticle whose moment can be represented by a Heisenberg 
spin with an easy axis anisotropy 0- To describe the 
dynamics of a Heisenberg spin, it is convenient to use 
the spherical coordinate system. The FPE in spherical 
coordinates 9 and Lp can be written in form as: 
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P — P{0, ip, t) is the probability density of the moment 
orientation. A and B are the so-called drift and diffusion 
coefficients respectively, defined as the ensemble mean of 
an infinitesimal change of 6 and Lp with respect to time 



The reduced stochastic LLG dynamical equation can 
be written as: 
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where m is the magnetic moment unit vector, a and 70 
are the damping and gyromagnetic constant respectively, 
h is the effective field normalized by the anisotropy field 
Hk = 2Ku/ HqMs, where is the anisotropy constant, 
/xq is the magnetic permeability and Ms is the saturation 
magnetization. The thermal field ht is introduced by 
Brown as a white noise term. The FPE corresponding 
to the LLG equation has been derived by Brown , and 
its factors are as follows: 




FIG. 1: Diagram of random walk step of length r and angle 
a to eg which define a spherical triangle ABC. 



rejection-free precession about an effective field. For the 
displacement about a cone, we pick a random vector ly- 
ing within a sphere of radius R to the original magnetic 
moment and then normalize the resulting vector. The 
precessional step vector, i.e. the displacement of mag- 
netic moment due to precession, is Am = — $ • m x h, 
where $ ^ 1 is a precesional step size to be determined. 
The probability q is chosen to be 1/2, which yields a 
near-optimal balance of efficiency and accuracy of our 
simulation. 

To calculate the FP coefficient A^''^ for the MC 
method, we obtain the ensemble mean of a small change 
of 9 in one Monte Carlo step. Contributions from random 
walk and precessional step are Ag^'^ = {A9)j.g^^^ /2 + 
A0prec/2, where the angle brackets denote the ensemble 
average. 

We first calculate the {A9)j.g^^^, where the angular dis- 
placement is defined by two random variables r and angle 
a, as shown in Fig. ^ After some geometrical analysis, 
we obtain: 



A 



LLG 



09 ~^ sm9 dip 



jiLLG 

B, 



LLG 



nLLG 
B, 

B, 



, I dE 
^ sm9 89 
2k' 
1 



+ k' cot ( 
1 dE 



sin^ 9 dip 



(3) 



LLG 



-2k' 



sm 




where in above equations, h' = ^ vM^{i+a'^) ' s' ^ h' /a, 
k' = h'/f], E is the total energy [^^|, V is the volume 
of the particle and /3 — {kBT)~^, ks is the Boltzmann 
constant and T is the temperature in Kelvin. 

We will now derive the FPE corresponding to our 
Monte Carlo method. For the MC method, we choose 
with probability q, to displace the magnetic moment 
within a small cone centered at the original magnetiza- 
tion direction, and with probability {1 — q) to perform a 
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Next, we require the probability for the displacement 
vector to be of size r (r < R) and angle a with re- 
spect to eg. This probability is given by Nowak et al. 
[la] as p{r) = 3VW~^/2ttR^. Based on the heat-bath 
Metropolis MC scheme, the acceptance rate is 
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where AE is the energy change in the random walk step. 
Thus, integrating over the projected surface of Fig.^ one 
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obtains {Ae)^.g^^^: 
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Next we calculate the other contribution from the pre- 
cessional step A^prec: 



A6'prec ^ eg ■ (-$ • m x h) 
1 $ 9E 
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In the above derivation, we have used the vector iden- 
tity a • (b X c) = (a X b) • c and h = — (Vm-E)/2i4r„F. 
Using Eqs. and Af^ becomes, 



Ar = |(cot.- 

(9) 

The other FPE factors can be obtained with the same 
procedure. 
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We can now compare the FPE factors corresponding to 
the Langevin (LLG) equation in Eq. Q , with those of the 
Metropolis MC method in Eqs. ^ and ifTUI) . Performing 
a term- wise comparison and omitting 0(R^) and higher 
order terms, we found that there is a one-to-one mapping 
between all FP terms of MC and LLG if: 
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Note that $ is in order of R? ^ thus we are justified in 
neglecting 0($^) terms in the above comparison between 
Eq. (PJ, and Eqs. © and From Eq. 1^1} , we obtain 
the time quantification factor of our hybrid Metropolis 
MC method, while Eq. H12() determines the precessional 
step size $. After taking into consideration the probabil- 
ity factor o, Eqs. (|ll|l and p2l) can be reduced to Nowak's 
results [l3 in the high damping case. 

To test the validity of Eqs. Hll|l and H12I) . we perform 
numerical calculations of the switching process for a mag- 
netic particle in which the easy axis is oriented at 7r/4 
to the applied field direction. All results are averaged 
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FIG. 2: Time dependence of magnetization along easy axis, 
for an isolated particle. KuV/kBT — 15, applied field h — 
0.42 tilted at 7r/4 relative to easy axis. Damping constant 
a = 0.5. 
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FIG. 3: Time dependence of magnetization along easy axis 
for interacting spin array. Periodic boundary conditions were 
used and KuV/ksT = 25, applied field h = 0.5 at a tilted 
angle of 7r/4 relative to the easy axis. Damping constant 
a = 1, exchange coupling strength J/Ku = 2 (Hamiltonian of 
an interacting system with exchange coupling strength J can 
be found, i.e. in Ref. 11^). R = 0.025 is used in the Monte 
Carlo simulation. Statistical error for the 10 x 10 lattice Monte 
Carlo simulation is shown in the inset. 



from a few thousand simulations. We consider the time 
evolution behavior of the mean magnetization compo- 
nent along the easy axis, and found a close convergence 
between our time-quantified MC method and the LLG 
equation (Fig.jSJ. In these calculations, we use R — 0.03 
for MC, and At = 0.001 for the LLG integration. We also 
apply our analytic results to 10 x 10, 20 x 20 and 40 x 40 
interacting spin array systems. For these simulations, 
R — 0.025 is used. We remark that a smaller step size 
R reduces statistical errors. We obtain a very good con- 
vergence between the MC and LLG results for all three 
arrays (Fig. O, especially so for the larger arrays. We 
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FIG. 4: Switching time v.s. damping constant a. 
KuV/ksT = 15, applied field h = 0.42 at a tilted angle of 
7r/4 relative to easy axis. Errorbars are smaller than the size 
of the symbols. Note that Nowak's method diverges from 
LLG equation at a < 2. 

believe this is due to the effects of self averaging. 

Note that our derivation of the FPE factors is appli- 
cable in a very general case. For instance, we do not 
require assumption of the system being in the vicinity of 
an energy minimum [l^ . The derivation also provides 
additional information, e.g., it explains mathematic ally 
why the Metropolis MC random walk method of Ref. [IJ 
fails to include the energy conservative precessional mo- 
tion. The FPE expression for the pure Metropolis MC 
does not contain terms corresponding to the g'-factor 



related terms of the LLG method [Eq. Q], which are 
precisely the terms which reflect the precessional part 
of the magnetization dynamics [T9l |. Thus, as shown in 
Fig. 21 we have successfully implemented the represen- 
tation of precessional motion in our MC method. We 
investigate the influence of damping constant on switch- 
ing time, where the switching time is defined as the time 
required for the magnetic moment to reach zero from the 
initial state. The precessional step size ^> guarantees a 
precise description of switching process even in the case of 
very low damping constant a, in which precessional mo- 
tion dominates the reversal process By contrast, the 
results obtained from the pure Metropolis MC method of 
Nowak et al. diverges significantly from that of the LLG 
equation at low a. 

To summarize, we have proposed a general method us- 
ing FPE to map MC to LLG and vice versa. We de- 
rived the drift and diffusion FP terms, corresponding to 
the MC method and compare them to the FP terms ob- 
tained from the LLG dynamical equation. By matching 
the terms in the drift and diffusion coefficients, we obtain 
a time quantification factor and the required precessional 
step size. The idea of using the FPE to link two stochas- 
tic methods (MC and Langevin methods) is general, and 
may be applied to other areas such as molecular dynam- 
ics in our future work. 
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